Metropolis Importance Sampling for Rugged Dynamical Variables 
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A funnel transformation is introduced, which acts recursively from higher towards lower temper- 
atures. It biases the a-priori probabilities of a canonical or generalized ensemble Metropolis sim- 
ulation, so that they zoom in on the global energy minimum, if a funnel exists indeed. A first, 
crude approximation to the full transformation, called rugged Metropolis one (RMi), is tested for 
Met-Enkephalin. At 300 K the computational gain is a factor of two and, due to its simplicity, RMi 
is well suited to replace the conventional Metropolis updating for these kind of systems. 
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To explain important aspects of protein folding, Bryn- 
gelson and Wolynes introduced a funnel picture [1], which 
is supported by various numerical results [2-4]. Never- 
theless, the understanding as well as the practical rele- 
vance of the funnel concept has remained somewhat lim- 
ited. The reason is that the funnel lives in the high- 
dimensional configuration space, while numerical studies 
have been confined to projections onto so called reaction 
coordinates, and there is no generic definition of a good 
reaction coordinate. Here I follow a different path and in- 
troduce a general funnel description from higher towards 
lower temperatures. It yields a powerful new method for 
designing Metropolis [5] weights. 

In protein models the energy E is a function of a num- 
ber of dynamical variables Vi, i = 1, . . . , n, whose fluc- 
tuations in the Gibbs canonical ensemble are described 
by a probability density (pd) . . . , w„; T), where T 

is the temperature. To be definite, we use in the follow- 
ing the all-atom energy function [6] ECEPP/2 (Empirical 
Conformational Energy Program for Peptides). Our dy- 
namical variables Vi are the dihedral angles, each chosen 
to be in the range — tt < Vi < tt, and the volume of the 
configuration space is K = (27r)". 

Let us define the support of a pd of the dihedral angles. 
Loosely speaking, the support of a pd is the region of con- 
figuration space where the protein wants to be. Mathe- 
matically, we define Rp to be the smallest sub- volume of 
the configuration space for which 
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holds. Here < p < 1 is a probability, which ought to 
be chosen close to one, e.g., p = 0.95. The free energy 
landscape at temperature T is called rugged, if the sup- 
port of the pd consists of many disconnected parts (this 
depends of course a bit on the adapted values for p and 
"many"). That a protein folds at room temperature, say 
300 K, into a unique native structure Vi, . . . means 



tion around this structure. We are now ready to formu- 
late the funnel picture in terms of pds. Let us choose a 
protein and consider for it a sequence of pds 

Pr{vi,...,Vn) = p{vi,...,Vn;Tr), r = 1,...,S , (2) 

which is ordered by the temperatures T^, namely 

Ti > T2 > ... > Tf . (3) 

The sequence (2) constitutes a protein funnel when, for 
a reasonable choice of the probability p and the temper- 
atures (3), the following holds: 

1. The pds are rugged. 

2. The support of a pd at lower temperature is con- 
tained in the support of a pd at higher temperature 



e.g. for p = 0.95, Ti ^ AOOK and T/ = 300 K. 



(4) 



that its pd p{vi, . 



i;300K) describes small fluctua- 



3. With decreasing temperatures the support Rp 
shrinks towards small fluctuations around the na- 
tive structure. 

Properties 2 and 3 are fulfilled for many systems of 
statistical physics, when some groundstate stands in for 
the native structure. The remarkable point is that they 
may still hold for certain complex systems with a rugged 
free energy landscape, i.e., with property 1 added. In 
such systems one finds typically local free energy min- 
ima, which are of negligible statistical importance at low 
temperatures, while populated at higher temperatures. 
In simulations at low temperature the problem of the 
molecular dynamics (for a review see [7]) as well as of the 
Metropolis [5] canonical ensemble approach is that the 
updating tends to get stuck in those local minima. On 
realistic simulation time scales this prevents convergence 
towards the native structure. On the other hand, the sim- 
ulations move quite freely at higher temperatures, where 
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the native structure is of negligible statistical weight. 
Nevertheless, the support of a protein pd may already 
be severely restricted, as we shall illustrate. The idea is 
to use a relatively easily calculable pd at a higher tem- 
perature to improve the performance of the simulation at 
a lower temperature. In the following we investigate this 
idea for the Metropolis algorithm. 

The Metropolis importance sampling would be per- 
fected, if we could propose new configurations {v'^} with 
their canonical pd p{v[, . . . ,v!^;T). Due to the funnel 
property 2 we expect that an estimate p{vi_, T') 
from some sufficiently close-by higher temperature T' > 
T will feed useful information into the simulation at tem- 
perature T. The potential for computational gains is 
large because of the funnel property 3. The suggested 
scheme for the Metropolis updating at temperature Tr 
is to propose new configurations {v^} with the pd (2) 
Pr-i{v[, . . . , v'^) and to accept them with the probability 
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This equation biases the a-priori probability of each di- 
hedral angle with an estimate of its pd from a higher 
temperature. In previous literature [8,9] such a biased 
updating has been used for the theory, where it is 
efficient to propose (f){i) at each lattice size i with its 
single-site probability. 

For our temperatures the ordering (3) is assumed. 
With the definition 'Pq{vi,. . . ,Vn) = (27r)~" the simula- 
tion at the highest temperature, Ti, is performed with 
the usual Metropolis algorithm. We have thus a recur- 
sive scheme, called rugged Metropolis (RM) in the fol- 
lowing. When p^_i{vi, . . . ,Vn) is always a usefull ap- 
proximation of pr{vi, . . . ,Vn), the scheme zooms in on 
the native structure, because the pd at Tf governs its 
fluctuations. 

To get things working, we need to construct an estima- 
tor . . . , f„; Tr) from the numerical data of the RM 
simulation at temperature T^. Although this is neither 
simple nor straightforward, a variety of approaches offer 
themselves to define and refine the desired estimators. In 
the following we work with the approximation 
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where the pl{vi, . . . ,Vn;Tr) are estimators of reduced 
one-variable pds defined by 



J — TT - / ■ 



p{vi,---,v„;T) . 
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The resulting algorithm, called RMi, appears to consti- 
tute the simplest RM scheme possible. Its implementa- 
tion is straightforward, as estimators of the one- variable 
reduced pds are easily obtained from the time series of 
a simulation. The computer time consumption of RMi 



is practically identical with the one of the conventional 
Metropolis algorithm. In the following RMi is used to 
demonstrate the correctness of our basic assumptions. 

The scope of this paper limits us to one illustra- 
tion. We rely on a simulation of the brain peptide Met- 
Enkephalin, because it is a numerically well studied sys- 
tem, which allows for comparison with results of the liter- 
ature [11-14]. Our Metropolis simulations are performed 
with a variant of SMMP [15] (Simple Molecular Mechan- 
ics for Proteins). We keep the u torsion angles uncon- 
strained and thus have 24 fully variable dihedral angles. 
In the previous literature the omega angles were either 
fixed to TT or restricted to [— 7r-|-7r/9,7r — 7r/9]. This leads 
to statistically significant differences in the energy, while 
the structural differences remain negligible. The perfor- 
mance of the RMi updating was tested at 300 K using 
input from a simulation at 400 K. The value of 300 K is 
chosen, because it is in the temperature range on which 
simulations of biological molecules eventually have to fo- 
cus. The temperature of 400 K is high enough so that 
the Metropolis algorithm is efficient as the autocorrela- 
tion times are small, while it is low enough to provide 
useful input for the 300 K simulation. 

At each simulation temperature a time series of 2^^ = 
131, 072 configurations is kept, in which subsequent con- 
figurations are separated by 32 sweeps. A sweep is de- 
fined by updating each dihedral angle once. Before start- 
ing with the measurements 2^^ = 262, 144 sweeps are per- 
formed for reaching equilibrium. Thus, the entire simula- 
tion at one temperature relies on 2^^ -|- 2^^ = 4, 456, 448 
sweeps. On a modern PC (1.9 GHz Athlon) this takes 
under 12 hours for the vacuum system and less than two 
days with the inclusion [16] of solvent effects. For each 
dihedral angle the acceptance rate of the Metropolis al- 
gorithm was monitored at run time and the integrated 
autocorrelation time Tint (see [10] for its definition) was 
calculated from the recorded time series. Values around 
0.5 are desirable for the acceptance rate, but the decisive 
quantity for the performance of an algorithm is the in- 
tegrated autocorrelation time. To achieve a pre-defined 
accuracy, the computer time needed is directly propor- 
tional to Tint - In the following results from vacuum sim- 
ulations are summarized. Computations which include 
solvent effects will be reported elsewhere [17]. 

In table I results for the energy and two dihedral angles 
are presented. Error bars are given in parenthesis. The 
acceptance rates are accurate to ±1 in their last digit. 
Using the SMMP [15] conventions, which differ from pre- 
vious literature [11-13], the angles are Gly-2 uj {vq) and 
Gly-3 cf) (wio). They are well-suited to illustrate impor- 
tant features of our approach. 

The Gly-3 (f) angle and four more angles (Gly-2 0, Gly- 
2 ijj, Phe-4 (f) and Phe-4 ip) exhibit very large autocor- 
relation times. In figure 1 estimates of the one- variable 
pds (6) for the Gly-3 <}) angle at 400 and 300 are 
depicted. The rugged nature of the distribution is al- 
ready obvious at 400 While the shapes of the other 
23 dihedral angle pds vary greatly, featuring from one to 
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TABLE I. Acceptance rates and integrated autocorrela- 
tions times for the energy E and, in the SMMP notation, 
the dihedral angles Gly-2 u) (ve) and Gly-3 4> (vio). 
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FIG. 1. Vacuum probability densities at 400 if and 300 AT 
for the Met-Enkephalin dihedral angle Gly-3 0. 

three local maxima, a ruggcdncss and our funnel prop- 
erty 2 are found for each case. Towards low temperatures 
a shrinking of the one- variable pds to their global energy 
minimum (GEM) implies property 3. The arrow in fig- 
ure 1 indicates the GEM value of this particular angle. 
Note that for the Metropolis algorithm the » regions 
do not constitute free energy barriers, because they can 
be jumped by a single Metropolis updating step. This is 
different for molecular dynamics simulations. 

Despite the similarities of the GIy-3 pds at 400 K 
and 300 K, there is a big increase of the Gly-3 4) inte- 
grated autocorrelation time. In the conventional, canon- 
ical Metropolis simulation it is by more than a factor of 
twenty, from 7.5 to 167. The RMi updating reduces this 
by a factor of 2.1, from 167 to 81. The integrated au- 
tocorrelations times of table I are given in units of 32 
sweeps, as this is the step-size of our time series data 
recorded. 



The acceptance rates for moves of the Gly-3 (j) angle 
are very small. As the support of its pds in figure 1 cov- 
ers more than 50% of the full range, this has to be due 
to correlations with other dihedral angles. To exhibit a 
rather distinct case of an angle with a low acceptance 
rate, results for one of the six oj angles are also included 
in table 1. For the conventional Metropolis simulation 
this angle has the same acceptance rate as the Gly-3 (j) 
angle (incidentally to all digits given). A look at the 
Gly-2 w pds reveals an obvious reason: They are nar- 
rowly peaked around the (identified) values ±7r, which is 
explained by the specific electronic hybridization of the 
CO-N peptide bond. Autocorrelations times arc about 
eight times smaller than for the Gly-3 (j) angle. Apply- 
ing RMi updating to the Gly-2 u pds cures entirely the 
problem of its low acceptance rate. At 300 K the increase 
is from 0.034 to 0.416 and similar numbers arc found for 
the other w angles. In contrast to that the Gly-3 ac- 
ceptance rate increases only to the modest value of 0.07, 
while the improvements of Tint arc kind of similar. For 
Gly-2 uj it is by a factor of 2.7 from 21 to 7.9. 

It is straightforward to combine the RMi updating 
with generalized ensemble methods (see [18] for a review). 
They are enabling techniques for studying equilibrium 
physics of complex systems at very low temperatures [19]. 
In the parallel tempering (PT) approach [20,21] two or 
more replica are simulated in parallel at distinct temper- 
atures and Metropolis exchanges of the temperatures are 
offered occasionally. Autocorrelations are reduced, when 
suitable excursions to higher temperatures become feasi- 
ble. In table I results for a PT simulation with replica 
at 400 K and 300 K are included. To compare with RMi, 
we use the integrated autocorrelation time of the energy, 
which characterizes the over-all performance. The PT 
method decreases Tint by a divisor of 2.5, from 50 to 20. 
When parallel nodes with a reasonably fast communica- 
tion are available, this is the gain in real time. For RMi 
the improvement factor is 1.9, from 50 to 26. Thus, PT 
outperforms RMi in real time, while RMi wins in the to- 
tal CPU time. Most remarkable is that the improvement 
factors multiply. Running PT with the RMi updating 
yields another factor of two. Tint is reduced from 20 to 10 
and, altogether, from 50 to 10. Note that the acceptance 
rates are not changed by PT. 

In previous simulations [13,3] it has been shown that 
considerably lower temperatures than 300 K arc needed 
to reduce the fluctuation of Met-Enkephalin in vacuum 
to fluctuation around its GEM. Here, I like to point out 
that each of our simulations at 300 K reaches the valley 
of attraction of the GEM suflScicntly often, so that the 
GEM can be found by local minimization. To be specific, 
the energy spectrum of figure 2 is obtained from our RMi 
simulation at 300 K in the following way: First wc isolate 
all configurations of the time series which are minima in 
the lower 10% q-tile of the energy distribution and sep- 
arated by an excursion of the time scries into the upper 
10% q-tile. In this way we obtain 1357 configurations and 
the SMMP minimizer is run on each of them. The prob- 
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FIG. 2. Spectrum of energy-minimized Met-Enkephalin 
configurations from the 300 K RMi simulation in vacuum. 



abilities of the resulting spectrum levels are plotted in 

figure 2. The GEM occurs at = -12.91 Kcal/mol with 
about 6% probability (107 times) followed by a conforma- 
tion at E = —10.92 Kcal/mol with about 4% probability. 
The frequency of finding the GEM for the other simula- 
tions of table I is approximately proportional to 7^"^ of 
the energy. 

Rugged distributions of the dynamical variables are 
typical for Metropolis simulations of proteins and RMi 
improves the importance sampling. It ought to become 
standard, as it yields a relevant gain in computer time for 
little extra efforts by the programmer. Met-Enkephalin 
is essentially solved by a simulation at 300 K. For some 
of its dihedral angles the RMi updating overcomes the 
problem of low acceptance rates entirely. For others the 
improvement remains more modest, because their low ac- 
ceptance rates are due to correlations with other angles. 
Only multi-variable moves can achieve importance sam- 
pling in such a situation. For the usual Metropolis algo- 
rithm the acceptance rates for multi-variable moves are 
practically zero. Here, we gained novel physical insight 
into the funnel picture, which provides us with a recipe 
on how to design multi-variable moves, so that the ac- 
ceptance rate is expected to stay reasonably large. One 
will start with the angles with the worst autocorrelations 
and construct their two-angle moves according to their 
reduced two-variable pds, the RM2 algorithm. Next, one 
may use three variables, and so on. In this way the out- 
come of our investigation of computational consequences 
of the funnel picture is that we do not expect a single, 
generically good Metropolis algorithm for protein simu- 
lations. Instead, we have developed a strategy to design 
an algorithm for each particular protein. This aspect of 
the RM approach promises advances for simulations of 
larger peptides. 
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